Multiple Numerical Assembly

Sections

Abstract

In this document, I want to provide a write-up of how this model is similar to and different from the versions previously implemented. This model is meant to provide assembly of multiple sites at the same time, that may, or may not, be connected in some fashion. Along the way, I will mention some of the input and output format that I am expect as documentation. At the end, I will be presenting some preliminary results. I especially want to use this as a vehicle to think about how to analyse those results.

Functions

library(dplyr)     # Data Manipulation

Warning: package ‘dplyr’ was built under R version 4.1.1

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal,
    union

library(ggplot2)   # 2-D Plot

Warning: package ‘ggplot2’ was built under R version 4.1.1

library(plotly)    # 3-D Plot

Warning: package ‘plotly’ was built under R version 4.1.1

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(ggfortify) # used for biplots of PCAs

Warning: package ‘ggfortify’ was built under R version 4.1.1

library(vegan)     # Ecological analysis mega-package

Warning: package ‘vegan’ was built under R version 4.1.1
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.1
Loading required package: lattice
This is vegan 2.5-7

library(RMTRCode2)

# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
    obj <- deparse(substitute(obj))
    if(exists(obj, envir = env)) {
        rm(list = obj, envir = env)
    }
}

Results, No Spatial Structure

First, we load up the preliminary results. In this case, we are loading a system in which 34 basal species and 66 consumer species form the pool for 10 unconnected environments. The pool and interaction matrix were assembled with the default parameters from Law and Morton’s 1996 work.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt-Result-Env10-None.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)

Events

In total, 9320 events were used in these environments, with the species and environment invasion both randomly assigned. The number of arrival and extinction events were controlled to both be half of this number. We chose this number due to the coupon collecting problem. In particular, we use the result that the probability of encountering each species is bounded: \[\text{Pr}(\text{Draws} < n \log_{e} n + c n) \rightarrow \exp(-\exp(-c)) \text{ as } n \rightarrow \infty\] where \(n\) is the number of species and \(c\) is a constant. For our purposes, we choose \(c = 5\) so that we have a probability of about \(99.3\%\) of seeing each species in each environment. In practice, we failed to observe 2 species-environment combinations. Notably, nearly every species had at least one successful invasion; 5 did not.

The initial abundance was set to be 4000 times the elimination threshold, in line with work on minimum viable populations (Traill et al. 2007). The elimination threshold is admittedly more arbitrary, since it sets an effective individual-area relationship. For this calculation, I set it to \(10^{-4}\), in line with our previous calculations. This is large enough to avoid numerical difficulties from precision, while being low enough to represent a decent sized region.

Since each population is assembling simultaneously, I chose to use exponential waiting times for the inter-species arrival and extinction times. Note that these rates are shared between species and environments, but arrival and extinction are fully independent of each other. Species and environment affected in each event was chosen uniformly at random. The question then is how to set the rate.

To set the rate in this case, I chose to set it to the largest eigenvalue magnitude of the per-environment interaction matrices. This magnitude corresponds to the strongest response we can see from the interaction matrix and, if the interaction matrix is a good approximation for the Jacobian around a stable fixed point (which is not guaranteed), indicates the characteristic time scale of the decay to equilibrium. Hence, (overall) arrival rates and (overall) extinction rates should happen on the same timescale as the (largest) dynamics in the system. Since there are 10 environments, we should then expect that 10 characteristic time scales, on average, should occur in between arrival events in the same environment.

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1, 
                                                 to = nrow(result$Abundance), 
                                                 by = 10), c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis

Every vertical line is a species introduction or extinction by neutral dynamics.

Alpha Diversity plots

Intro

Perhaps more intriguing might be some sense of the biodiversity that we have in each system. We break the abundance results up by environment, then calculate the number of non-zero abundance curves at each time point. We also calculate the Shannon entropy (reminder: higher entropy means more uncertainty which means a flatter distribution).

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Mean = mean(Richness),
      SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
        Species, collapse = ", "), split = ", ", fixed = TRUE))))),
      Gamma = unlist(lapply(strsplit(
        SpeciesTotal, split = ", ", fixed = TRUE), length))
    ) %>% tidyr::pivot_longer(
      cols = c(Mean, Gamma), 
      names_to = "Aggregation",
      values_to = "Richness"
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness,
    linetype = Aggregation
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

So richness hovers around a similar regime throughout the majority of the simulation. Note in this plot that we have emphasised one environmental curve and superimposed the mean in black. We do manage to reach heights of 11 species in one environment, but these heights are shortlived. Instead, we seem to observe a (time and environment averaged) value:

mean(Diversity$Richness)

[1] 4.458893

(If we consider the first 10,000 time units as burn-in, we instead see a value of

mean(Diversity$Richness[Diversity$Time > 10000])

[1] 4.587217

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 8043 row(s) containing missing values (geom_path).
Warning: Removed 2297 row(s) containing missing values (geom_path).

Entropy has a similar, but highly erratic, behaviour. If one tries to follow one of the entropy curves, then one sees that they have fairly substantial periods of almost smooth behaviour followed by suddenly very noisy behaviour, and noise seems to be the dominant mode if one tries to examine the low alpha environment in the background. There are some easy to make predictions. Extinctions reduce the entropy in the system, as you become more certain about what remains. Analogously, arrivals increase the entropy. We can probably better see the relationship beyond these principles by plotting entropy against richness and connecting observations by environment and time.

Entropy-Richness

# ggplot2::ggplot(
#   Diversity %>% dplyr::filter(Environment < 4), 
#   ggplot2::aes(
#     x = Richness,
#     y = Entropy,
#     group = Environment,
#     color = Time
#   )
# ) + ggplot2::geom_path(
# ) + ggplot2::guides(
#   alpha = "none"
# )

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

It seems quite hard to tell, but there does not appear to be any particular orientation (clockwise, counter-clockwise) or similar pattern here.

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 10986 row(s) containing missing values (geom_path).
Warning: Removed 2831 row(s) containing missing values (geom_path).

Evenness helps highlight that this is a high variance process but with a relatively constrained mean.

Environment Diversity

We can, of course, flip the idea on its head. Instead of examining the diversity of species within environments, we can look at the diversity of environments occupied by species. Since very few species end up occupying environments, we just look at richness. Unfortunately, this is quite a memory exhaustive task.

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

One immediately interesting trend here is that very few species are present across more than 5 environments at a given time. Indeed, only

unique(EnvDiversity$Species[EnvDiversity$Richness > 5])

[1] 12 14 22 30

are ever present in more than 5 environments at once. We can also examine how long these periods occur for by species by tabulation. We block times so that entries are all of the same unit length, and round the average richness during the given time unit.

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species     0     1     2     3     4     5     6     7
    1   65970 14384   439     0     0     0     0     0
    2   34776 31378 11352  3287     0     0     0     0
    3   63206 13000  1767  2820     0     0     0     0
    4   46367 21078 11443  1617   288     0     0     0
    5   30648 31553 16542  1341   709     0     0     0
    6   72061  8732     0     0     0     0     0     0
    7   75522  5271     0     0     0     0     0     0
    8   64759 15508   526     0     0     0     0     0
    9   77919  2874     0     0     0     0     0     0
    10  17817 38679 20192  4105     0     0     0     0
    11  52666 24892  3235     0     0     0     0     0
    12   9669  2676  9781 26641 18597  9922  2294  1213
    13  64640 15703   450     0     0     0     0     0
    14  12094 21314 14556 16445 15148   931   305     0
    15  68564 10732  1497     0     0     0     0     0
    16  69690 10903   200     0     0     0     0     0
    17  30636 33714 15505   938     0     0     0     0
    18  69214 11322   257     0     0     0     0     0
    19  62460 17599   734     0     0     0     0     0
    20  74867  3685  1849   392     0     0     0     0
    21  47906 25702  4185  3000     0     0     0     0
    22   3180 13477 20151 21766 16740  5305   174     0
    23  65137 14940   716     0     0     0     0     0
    24  65792 12313  2688     0     0     0     0     0
    25  66085 11109  3599     0     0     0     0     0
    26  41654 19402 19201   536     0     0     0     0
    27  48808 20447 10279  1259     0     0     0     0
    28  68983 11810     0     0     0     0     0     0
    29  34631 26993 17044  2125     0     0     0     0
    30   5561 19065 20682 11364 11597  7963  3301  1260
    31  65019 15665   109     0     0     0     0     0
    32  19327 26015 29600  5851     0     0     0     0
    33  15128 23405 24345 10059  7856     0     0     0
    34  76432  4361     0     0     0     0     0     0
    35  73829  6080   884     0     0     0     0     0
    36  71050  9743     0     0     0     0     0     0
    37  38199 22369 15663  4486    76     0     0     0
    38  28752 35327 16714     0     0     0     0     0
    39  77008  3785     0     0     0     0     0     0
    40  49953 27129  3711     0     0     0     0     0
    41  79838   955     0     0     0     0     0     0
    42  78312  2481     0     0     0     0     0     0
    43  32813 23892 17056  6661   371     0     0     0
    44  39700 23980 15554  1559     0     0     0     0
    45  80793     0     0     0     0     0     0     0
    46  48461 27669  4663     0     0     0     0     0
    47  12053 36783 14826 15182  1949     0     0     0
    48  63151 16980   662     0     0     0     0     0
    49  73125  7668     0     0     0     0     0     0
    50  79860   933     0     0     0     0     0     0
    51  80255   538     0     0     0     0     0     0
    52  73946  6847     0     0     0     0     0     0
    53  77581  3212     0     0     0     0     0     0
    54  76399  4394     0     0     0     0     0     0
    55  55217 23951  1625     0     0     0     0     0
    56  72564  8229     0     0     0     0     0     0
    57  57443 22644   706     0     0     0     0     0
    58  78988  1805     0     0     0     0     0     0
    59  35315 37942  7536     0     0     0     0     0
    60  67881 12912     0     0     0     0     0     0
    61  52226 24749  3818     0     0     0     0     0
    62  62517 16663  1613     0     0     0     0     0
    63  46937 31441  2415     0     0     0     0     0
    64  68974 11213   606     0     0     0     0     0
    65  80793     0     0     0     0     0     0     0
    66  67130 12561  1102     0     0     0     0     0
    67  37887 27036 11834  4036     0     0     0     0
    68  51139 27047  2205   402     0     0     0     0
    69  61791 16492  2510     0     0     0     0     0
    70  73231  7354   208     0     0     0     0     0
    71  22117 22512 20026 14172  1966     0     0     0
    72  10060 26081 33450  7590  3612     0     0     0
    73  73696  7097     0     0     0     0     0     0
    74  80793     0     0     0     0     0     0     0
    75  77528  3008   257     0     0     0     0     0
    76  34905 34766 10606   516     0     0     0     0
    77  54799 21256  4116   622     0     0     0     0
    78  57557 21260  1976     0     0     0     0     0
    79  46966 28430  5397     0     0     0     0     0
    80  46394 30478  3921     0     0     0     0     0
    81  80793     0     0     0     0     0     0     0
    82  72153  8640     0     0     0     0     0     0
    83  57732 21901  1160     0     0     0     0     0
    84  73061  7732     0     0     0     0     0     0
    85  78108  2685     0     0     0     0     0     0
    86  64841 15952     0     0     0     0     0     0
    87  75894  4488   411     0     0     0     0     0
    88  77088  3705     0     0     0     0     0     0
    89  48328 32465     0     0     0     0     0     0
    90  80793     0     0     0     0     0     0     0
    91  36647 30720 11100  1868   458     0     0     0
    92  54513 24626  1654     0     0     0     0     0
    93  73946  6847     0     0     0     0     0     0
    94  69274 10406  1113     0     0     0     0     0
    95  43213 20293 14643  2644     0     0     0     0
    96  65015 15778     0     0     0     0     0     0
    97  75428  5365     0     0     0     0     0     0
    98  65968 13004  1821     0     0     0     0     0
    99  78629  2164     0     0     0     0     0     0
    100 52405 15982  9604  2802     0     0     0     0

Beta Diversity

Principal Components Analysis

(For the desktop, the amount of data we generated is too much, so we need to reduce the amount of data we use. To do so we average over time blocks, here of length 100.)

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

We can perform a PCA and see if are data can be summarised by a small number of dimensions. As we shall see, constraining to the first 25 principal components does not harm the system much.

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

There is not actually a lot of dependence within the system it appears. Consider the amount of variance explained by the first six principal components (ordered, as is tradition, by amount of variation explained).

head(summary(PCA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.03749 0.03273 0.03095 0.02751 0.02501 0.02482 

The amount of explained variation is as follows.

sum(summary(PCA)$importance[2, ])

[1] 0.99994

The traditional biplot follows, but despite the seeming presence of patterns, so little of the variance is explained that is probably not worth further examination.

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

We can try again with presence-absence data instead.

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.05537 0.04584 0.04483 0.03901 0.03578 0.03410 

So it is a bit better, but not by much.

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

I should note that this is not a failure of the method persay. What this says to me is that dimension reduction of the system cannot reduce the system to a human-readable set of descriptors. The system is still being reduced (from 100 Species * 10 Environments to 25 Components to cover

sum(summary(PCA)$importance[2, ])

[1] 0.99994

of the variation). My initial hypothesis for when the systems are coupled is that, as coupling increases, the number of principal components should decrease, since one system’s changes will be better predictors of the next system’s changes. (An interesting related question; do the number of principal components correlate with the amount of biodiversity in the system?

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Since trying to reduce the system dimensionally does not work (well enough) a next attempt might be to consider how the pair-wise beta diversity changes over the course of the simulation. Initial problems include that there are quite a few measures of beta diversity as well as that the (square of the) number of environments will determine how many entries we need to consider.

To try this, we are going to reorganise our data into data frame with attributes of environment, time, and traits/species.

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

One perhaps interesting idea is to try to group the environments by cluster by considering the abundances (or presence-absence) of the species present as traits. The question then is whether the environments have a tendency to attract to specific points or if they wander around each other without relation. (The latter is more neutral.) If they appear to be convergent (which so far would seem to disagree mostly with the alpha diversity analyses), then that would imply that dynamics determine the majority of the system’s fate, while if they instead wander more randomly, then they would appear to be dominated by the neutral mechanisms. (Of course, this is affected by the rate of the neutral mechanisms, so it exists on a definite continuum.) (Note, of course, that cluster analysis usually includes something like PCA, so we would not necessarily expect a different result.)

My working hypothesis for how to use the distance measures is that we need a metric measure (Kindt and Coe, 2005) with any identification scrubbed off (site, time collected). Kindt and Coe suggest that taking the square root of Bray-Curtis makes it metric and thus usable for mathematical distance analyses (such as hierarchical clustering). The vegan package recomends using Jaccard instead of Bray-Curtis, for the same reason.

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

Note that are major breakdowns extremely early yielding a large number of clusters but there are expected to be 10 environments (if historical contingency mattered most) or a small number of clusters (if they are converging). The large number of clusters says to me that the neutral dynamics are jumping rapidly from cluster to cluster.

We can try again with presence absence instead.

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)

Results, Line Spatial Structure

ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
ifrm(EnvironmentDistancePA)

We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1, 
                                                 to = nrow(result$Abundance), 
                                                 by = 10), c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis

Every vertical line is a species introduction or extinction by neutral dynamics. Note that, due to dispersal, these should lose some of their abundance to spreading, but if that spreading is faster than the population can recoup its losses, than it can invade and be wiped out (since invadability does not consider spatial dynamics, only local ecological dynamics). The primary notes then are that there are a large number more events appearing to occur (due to the spatial arrangement), but the system as a whole is stabler (since the species that get knocked out can be recaptured from space as well).

Alpha Diversity plots

Intro

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Mean = mean(Richness),
      SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
        Species, collapse = ", "), split = ", ", fixed = TRUE))))),
      Gamma = unlist(lapply(strsplit(
        SpeciesTotal, split = ", ", fixed = TRUE), length))
    ) %>% tidyr::pivot_longer(
      cols = c(Mean, Gamma), 
      names_to = "Aggregation",
      values_to = "Richness"
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness,
    linetype = Aggregation
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1123 row(s) containing missing values (geom_path).
Warning: Removed 110 row(s) containing missing values (geom_path).

Entropy-Richness

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1440 row(s) containing missing values (geom_path).
Warning: Removed 136 row(s) containing missing values (geom_path).

Environment Diversity

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species     0     1     2     3     4     5     6     7     8     9    10
    1   80793     0     0     0     0     0     0     0     0     0     0
    2   80624     0     1     1     2     1     1     2     0     2   159
    3   80793     0     0     0     0     0     0     0     0     0     0
    4   80047     0     0     2     0     1     0     0     1     0   742
    5   79905     0     6     0     6     2     5     1     0    11   857
    6   80793     0     0     0     0     0     0     0     0     0     0
    7   80793     0     0     0     0     0     0     0     0     0     0
    8   80793     0     0     0     0     0     0     0     0     0     0
    9   80793     0     0     0     0     0     0     0     0     0     0
    10  80793     0     0     0     0     0     0     0     0     0     0
    11  80793     0     0     0     0     0     0     0     0     0     0
    12    680     0     0     0     0     0     0     0     0     1 80112
    13  80639     1     1     1     1     1     1     2     1     1   144
    14     84     0     0     0     0     0     0     0     0     0 80709
    15  80793     0     0     0     0     0     0     0     0     0     0
    16  80793     0     0     0     0     0     0     0     0     0     0
    17  80793     0     0     0     0     0     0     0     0     0     0
    18  80793     0     0     0     0     0     0     0     0     0     0
    19  80793     0     0     0     0     0     0     0     0     0     0
    20  80793     0     0     0     0     0     0     0     0     0     0
    21  80793     0     0     0     0     0     0     0     0     0     0
    22     62     0     0     0     0     0     0     0     0     1 80730
    23  80793     0     0     0     0     0     0     0     0     0     0
    24  80546     0     1     0     0     0     0     1     0     0   245
    25  80793     0     0     0     0     0     0     0     0     0     0
    26  80793     0     0     0     0     0     0     0     0     0     0
    27  80793     0     0     0     0     0     0     0     0     0     0
    28  80793     0     0     0     0     0     0     0     0     0     0
    29  80721     1     0     0     1     0     0     0     0     1    69
    30  69391     9     2     8     1     8     5     2     8    12 11347
    31  80793     0     0     0     0     0     0     0     0     0     0
    32    500     0     0     0     0     0     1     0     0     1 80291
    33  66491     5    15    10    14     2    14     5    14    17 14206
    34  80793     0     0     0     0     0     0     0     0     0     0
    35  80332     1     0     1     1     1     1     1     0     2   453
    36  80793     0     0     0     0     0     0     0     0     0     0
    37  46279     1     1     2     1     1     3     1     2     4 34498
    38  80793     0     0     0     0     0     0     0     0     0     0
    39  80793     0     0     0     0     0     0     0     0     0     0
    40  80793     0     0     0     0     0     0     0     0     0     0
    41  80793     0     0     0     0     0     0     0     0     0     0
    42  80793     0     0     0     0     0     0     0     0     0     0
    43  27630     2     2     0     4     0     0     1     4     7 53143
    44  60232     1     8     3     3     7     2     7     5    14 20511
    45  80793     0     0     0     0     0     0     0     0     0     0
    46  80793     0     0     0     0     0     0     0     0     0     0
    47  52191     5     8     3     7     3     7     5     9    12 28543
    48  80793     0     0     0     0     0     0     0     0     0     0
    49  80793     0     0     0     0     0     0     0     0     0     0
    50  80793     0     0     0     0     0     0     0     0     0     0
    51  80793     0     0     0     0     0     0     0     0     0     0
    52  80793     0     0     0     0     0     0     0     0     0     0
    53  80793     0     0     0     0     0     0     0     0     0     0
    54  80793     0     0     0     0     0     0     0     0     0     0
    55  78735     1     1     2     0     1     2     0     0     2  2049
    56  80793     0     0     0     0     0     0     0     0     0     0
    57  79020     0     0     1     0     0     0     0     1     0  1771
    58  80793     0     0     0     0     0     0     0     0     0     0
    59  80793     0     0     0     0     0     0     0     0     0     0
    60  80793     0     0     0     0     0     0     0     0     0     0
    61  80466     0     1     0     0     1     1     1     0     1   322
    62  80793     0     0     0     0     0     0     0     0     0     0
    63  80793     0     0     0     0     0     0     0     0     0     0
    64  80793     0     0     0     0     0     0     0     0     0     0
    65  80793     0     0     0     0     0     0     0     0     0     0
    66  80793     0     0     0     0     0     0     0     0     0     0
    67  80793     0     0     0     0     0     0     0     0     0     0
    68  79023     1     5     4     7     5     6     4     4     4  1730
    69  78135     2     2     1     1     1     6     1     2     5  2637
    70  80793     0     0     0     0     0     0     0     0     0     0
    71   1394     0     0     0     0     0     0     0     0     1 79398
    72    292     1     0     0     0     0     1     0     1     0 80498
    73  80793     0     0     0     0     0     0     0     0     0     0
    74  80793     0     0     0     0     0     0     0     0     0     0
    75  80793     0     0     0     0     0     0     0     0     0     0
    76  78676     1     4     1     4     0     6     1     3     6  2091
    77  79765     0     0     1     0     0     1     0     0     1  1025
    78  79253     0     0     1     0     1     0     0     0     0  1538
    79  80793     0     0     0     0     0     0     0     0     0     0
    80  80793     0     0     0     0     0     0     0     0     0     0
    81  80793     0     0     0     0     0     0     0     0     0     0
    82  80793     0     0     0     0     0     0     0     0     0     0
    83  80793     0     0     0     0     0     0     0     0     0     0
    84  80793     0     0     0     0     0     0     0     0     0     0
    85  80793     0     0     0     0     0     0     0     0     0     0
    86  64862     3     6     4     7     4     6     3     8     6 15884
    87  80793     0     0     0     0     0     0     0     0     0     0
    88  80793     0     0     0     0     0     0     0     0     0     0
    89  80793     0     0     0     0     0     0     0     0     0     0
    90  80793     0     0     0     0     0     0     0     0     0     0
 [ reached getOption("max.print") -- omitted 10 rows ]

Beta Diversity

Principal Components Analysis

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.19145 0.11089 0.10844 0.07943 0.05993 0.04530 

sum(summary(PCA)$importance[2, ])

[1] 0.99996

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.21750 0.12185 0.08092 0.05943 0.05215 0.04468 

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

sum(summary(PCA)$importance[2, ])

[1] 0.99996

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)

Results, Slow Line Spatial Structure

ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)

We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt1e+06-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[1:(length(toRemove)/2)][toRemove[1:(length(toRemove)/2)]] <- 0
result$Abundance[(length(toRemove)/2 + 1):length(toRemove)][toRemove[(length(toRemove)/2 + 1):length(toRemove)]] <- 0
ifrm(toRemove)

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1, 
                                                 to = nrow(result$Abundance), 
                                                 by = 10), c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis

Alpha Diversity plots

Intro

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Mean = mean(Richness),
      SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
        Species, collapse = ", "), split = ", ", fixed = TRUE))))),
      Gamma = unlist(lapply(strsplit(
        SpeciesTotal, split = ", ", fixed = TRUE), length))
    ) %>% tidyr::pivot_longer(
      cols = c(Mean, Gamma), 
      names_to = "Aggregation",
      values_to = "Richness"
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness,
    linetype = Aggregation
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1855 row(s) containing missing values (geom_path).
Warning: Removed 259 row(s) containing missing values (geom_path).

Entropy-Richness

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 2594 row(s) containing missing values (geom_path).
Warning: Removed 329 row(s) containing missing values (geom_path).

Environment Diversity

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species     0     1     2     3     4     5     6     7     8     9    10
    1   80793     0     0     0     0     0     0     0     0     0     0
    2   73294  1308   345  5846     0     0     0     0     0     0     0
    3   60300 12419   913  7161     0     0     0     0     0     0     0
    4   71970  7391   681   751     0     0     0     0     0     0     0
    5   77639  2117    30  1007     0     0     0     0     0     0     0
    6   80793     0     0     0     0     0     0     0     0     0     0
    7   80793     0     0     0     0     0     0     0     0     0     0
    8   80793     0     0     0     0     0     0     0     0     0     0
    9   80793     0     0     0     0     0     0     0     0     0     0
    10  62872  8776  1697  7448     0     0     0     0     0     0     0
    11  80505   246    42     0     0     0     0     0     0     0     0
    12    680    57     1    67    50    47    85    98   149   131 79428
    13  80793     0     0     0     0     0     0     0     0     0     0
    14     84    17   540   961  2714  3766  9410  3966  6741 12115 40479
    15  71934  2740    85  6034     0     0     0     0     0     0     0
    16  80793     0     0     0     0     0     0     0     0     0     0
    17  47555  7152   340  5033  5534  8193  1609  1541  3836     0     0
    18  80499   294     0     0     0     0     0     0     0     0     0
    19  80686   107     0     0     0     0     0     0     0     0     0
    20  80793     0     0     0     0     0     0     0     0     0     0
    21  79313  1480     0     0     0     0     0     0     0     0     0
    22     62    14     0    32     0    36     0    31     0  2888 77730
    23  80793     0     0     0     0     0     0     0     0     0     0
    24  80793     0     0     0     0     0     0     0     0     0     0
    25  80793     0     0     0     0     0     0     0     0     0     0
    26  80407   386     0     0     0     0     0     0     0     0     0
    27  80793     0     0     0     0     0     0     0     0     0     0
    28  80793     0     0     0     0     0     0     0     0     0     0
    29  79711   114   387   581     0     0     0     0     0     0     0
    30  23183 15781 13353 12433 16043     0     0     0     0     0     0
    31  80793     0     0     0     0     0     0     0     0     0     0
    32   4784 13308 14738  6522  8650 13321  6596  6339  2225  3634   676
    33  34246 12708 12031  8137 10278  2444   722   227     0     0     0
    34  80685   108     0     0     0     0     0     0     0     0     0
    35  80600   193     0     0     0     0     0     0     0     0     0
    36  56841 12740 11212     0     0     0     0     0     0     0     0
    37   7020 40060 26246  7467     0     0     0     0     0     0     0
    38  73593  7200     0     0     0     0     0     0     0     0     0
    39  80793     0     0     0     0     0     0     0     0     0     0
    40  78481  2312     0     0     0     0     0     0     0     0     0
    41  80271   522     0     0     0     0     0     0     0     0     0
    42  73553  7240     0     0     0     0     0     0     0     0     0
    43  18796 19884 26680 11558  3875     0     0     0     0     0     0
    44  30628 14552  5154 12331 12054  6074     0     0     0     0     0
    45  77510  3283     0     0     0     0     0     0     0     0     0
    46  70547  7937   160  2149     0     0     0     0     0     0     0
    47    979  1092 10511 17784 13851 20944 11456  4176     0     0     0
    48  52054 21980  6759     0     0     0     0     0     0     0     0
    49  70981  9812     0     0     0     0     0     0     0     0     0
    50  78655  2138     0     0     0     0     0     0     0     0     0
    51  79858   935     0     0     0     0     0     0     0     0     0
    52  69649  2937   894  5160  2153     0     0     0     0     0     0
    53  79424  1369     0     0     0     0     0     0     0     0     0
    54  78212  1554  1027     0     0     0     0     0     0     0     0
    55  65177 15616     0     0     0     0     0     0     0     0     0
    56  69587 11206     0     0     0     0     0     0     0     0     0
    57  52153 20987  3657  3988     8     0     0     0     0     0     0
    58  71255  9538     0     0     0     0     0     0     0     0     0
    59  56304 24489     0     0     0     0     0     0     0     0     0
    60  74416  6377     0     0     0     0     0     0     0     0     0
    61  67723 12110   960     0     0     0     0     0     0     0     0
    62  66212  4509 10072     0     0     0     0     0     0     0     0
    63  75952  4841     0     0     0     0     0     0     0     0     0
    64  76302  4491     0     0     0     0     0     0     0     0     0
    65  80600   193     0     0     0     0     0     0     0     0     0
    66  76287  4506     0     0     0     0     0     0     0     0     0
    67  51275 22463  7055     0     0     0     0     0     0     0     0
    68  26766 14846  7523 22838  7263  1557     0     0     0     0     0
    69  70929  8209  1600    55     0     0     0     0     0     0     0
    70  75050  5039   704     0     0     0     0     0     0     0     0
    71   7364  4049 25788 11967 18316 13309     0     0     0     0     0
    72    293  6015 31082 15185  8897  9826  6663  2832     0     0     0
    73  65220 15573     0     0     0     0     0     0     0     0     0
    74  76900  3893     0     0     0     0     0     0     0     0     0
    75  78673  2120     0     0     0     0     0     0     0     0     0
    76  60714 19313   662   104     0     0     0     0     0     0     0
    77  66756  3870  5716   138   704   173   375   351  1496  1214     0
    78  68931 11862     0     0     0     0     0     0     0     0     0
    79  74538  6255     0     0     0     0     0     0     0     0     0
    80  70951  9200   642     0     0     0     0     0     0     0     0
    81  79192  1601     0     0     0     0     0     0     0     0     0
    82  80444   349     0     0     0     0     0     0     0     0     0
    83  79662  1131     0     0     0     0     0     0     0     0     0
    84  73288  7505     0     0     0     0     0     0     0     0     0
    85  78008  2785     0     0     0     0     0     0     0     0     0
    86   8225 25572 17279 24066  5651     0     0     0     0     0     0
    87  79538  1255     0     0     0     0     0     0     0     0     0
    88  80538   255     0     0     0     0     0     0     0     0     0
    89  38907 34291  7595     0     0     0     0     0     0     0     0
    90  80793     0     0     0     0     0     0     0     0     0     0
 [ reached getOption("max.print") -- omitted 10 rows ]

Beta Diversity

Principal Components Analysis

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.06245 0.05166 0.04374 0.03644 0.03204 0.02918 

sum(summary(PCA)$importance[2, ])

[1] 1.00007

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.09404 0.07691 0.05775 0.05016 0.04514 0.03782 

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

sum(summary(PCA)$importance[2, ])

[1] 1.00007

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)